{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Mesh generation for tutorial 06, cases 1a, 2a, 3a, 4a-b, 5, 6, 8\n",
    "\n",
    "This file generates the mesh which is used in the following examples:\n",
    "* 1a_poisson\n",
    "* 2a_advection_diffusion_reaction\n",
    "* 3a_advection_diffusion_reaction_neumann_control\n",
    "* 4a_poisson_dirichlet_control\n",
    "* 4b_poisson_neumann_control_boundary_observation\n",
    "* 5_stokes\n",
    "* 6_navier_stokes\n",
    "* 8_navier_stokes_neumann_control"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from mpi4py import MPI\n",
    "from dolfinx import MeshTags, UnitSquareMesh\n",
    "from dolfinx.io import XDMFFile\n",
    "from dolfinx.mesh import locate_entities_boundary"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "assert MPI.COMM_WORLD.size == 1, \"This mesh generation notebook is supposed to be run in serial\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Create mesh"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mesh = UnitSquareMesh(MPI.COMM_WORLD, 32, 32)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Create (trivial) subdomains"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "num_cells = mesh.topology.index_map(mesh.topology.dim).size_local\n",
    "subdomains = MeshTags(mesh, mesh.topology.dim,\n",
    "                      np.arange(0, num_cells), np.zeros((num_cells,), dtype=np.intc))\n",
    "subdomains.name = \"subdomains\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Create boundaries"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def bottom(x):\n",
    "    return abs(x[1] - 0.) < np.finfo(float).eps\n",
    "\n",
    "\n",
    "def left(x):\n",
    "    return abs(x[0] - 0.) < np.finfo(float).eps\n",
    "\n",
    "\n",
    "def top(x):\n",
    "    return abs(x[1] - 1.) < np.finfo(float).eps\n",
    "\n",
    "\n",
    "def right(x):\n",
    "    return abs(x[0] - 1.) < np.finfo(float).eps\n",
    "\n",
    "\n",
    "bottom.id = 1\n",
    "left.id = 2\n",
    "top.id = 3\n",
    "right.id = 4\n",
    "\n",
    "boundaries_entities = dict()\n",
    "boundaries_values = dict()\n",
    "for boundary in (bottom, left, top, right):\n",
    "    boundaries_entities[boundary.id] = locate_entities_boundary(\n",
    "        mesh, mesh.topology.dim - 1, boundary)\n",
    "    boundaries_values[boundary.id] = np.full(\n",
    "        boundaries_entities[boundary.id].shape, boundary.id, dtype=np.intc)\n",
    "boundaries = MeshTags(mesh, mesh.topology.dim - 1,\n",
    "                      np.hstack(list(boundaries_entities.values())),\n",
    "                      np.hstack(list(boundaries_values.values())))\n",
    "boundaries.name = \"boundaries\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Save"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with XDMFFile(MPI.COMM_WORLD, \"square.xdmf\", \"w\") as output:\n",
    "    output.write_mesh(mesh)\n",
    "    output.write_meshtags(subdomains)\n",
    "    output.write_meshtags(boundaries)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython"
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
